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This paper presents an approach to shape an aircraft to equivalent area based ob- 
jectives using the discrete adjoint approach. Equivalent areas can be obtained either 
using reversed augmented Burgers equation or direct conversion of off-body pressures 
into equivalent area. Formal coupling with CFD allows computation of sensitivities of 
equivalent area objectives with respect to aircraft shape parameters. The exactness of 
the adjoint sensitivities is verified against derivatives obtained using the complex step 
approach. This methodology has the benefit of using designer-friendly equivalent areas 
in the shape design of low-boom aircraft. Shape optimization results with equivalent 
area cost functionals are discussed and further refined using ground loudness based 
objectives. 


Nomenclature 


Co Ambient speed of sound, m/s. 

C v Dimensionless dispersion. 

D Vector of design variables. 

G Ray tube area, m 2 . 

k n Scaling factor due to ray-tube speading and stratification. 

L Lagrangian. 

l n Cost function for adjoint calculation. 

Dispersion parameter. 

N Number of steps during propagation. 

p Pressure waveform during propagation. 

p t Target ground signature. 

q,r,t Intermediate pressure waveforms. 

t Retarded time. 

A n , B n Matrices during first relaxation process. 

A%, Blf Matrices during second relaxation process. 

Matrices during absorption process. 

P 1 + V- 

5 Diffusion parameter, 

r Dimensionless thermo-viscous parameter. 

* Senior Research Engineer, National Institute of Aerospace. 
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Ratio of specific heats, 1.4. 


A„, p n , 7o,n, 7i,n Adjoint vectors. 
ujq Angular frequency. 

Po Ambient density. 

Tu Dimensionless time for each relaxation mode. 

t' Intermediate retarded time coordinate. 

Q v Dimensionless relaxation time parameter. 

I. Introduction and Motivation 

Development of novel and useful methods for sonic boom mitigation of civil supersonic aircraft remains 
one of the most important steps in the conceptual and preliminary design. Since the 1960’s, researchers 
realized 1,2,3 the importance of aircraft shaping in reducing the sonic boom impact. The Shaped Sonic 
Boom Demonstrator (SSBD) 4 program verified, via flight testing, that aircraft shaping is an effective 
strategy for changing the boom signature on the ground. Ever since this achievement, there has been a 
renewed interest towards achieving better designs aimed at reducing the strength of the boom footprint. 

While there have been studies 5,6 showing promise towards achieving sonic boom mitigation without 
the use of sensitivity information, design approaches based on sensitivity of the chosen cost functional 
with respect to aircraft shape offer an effective and theoretically sound way to reduce the adverse impact 
of sonic boom. Adjoint-based methods offer an efficient way of computing sensitivities of various aerody- 
namic quantities to many shape design parameters. Several studies have demonstrated the capabilities 
of adjoint-based methods to optimize near-field sonic boom 7,8 or ground loudness 9 cost functionals. 

In the current study, aerodynamic analysis is performed using the FUN3D CFD code developed at 
NASA Langley. FUN3D provides discretely consistent adjoint capabilities for sensitivity analysis. It has 
been extensively used to perform adjoint-based mesh adaptation 10,11,12,13 and design optimization, 14,15 
including optimization of near-field sonic-boom waveforms. In the present paper, we are looking beyond 
the current state-of-the-art approaches and seek to optimize and match designer-friendly equivalent area 
distributions rather than a near-field target. This capability represents a third avenue to mitigate sonic 
boom, the other two being: a) near-field target matching 7,8 and b) using ground-based objectives. 16,9 
The reason for developing an adjoint for equivalent area matching can be summarized briefly as follows: 
analysis and design based on equivalent areas remains an attractive option to aircraft designers due to 
physics that are aligned with engineering intuition. This paper, however, does not attempt to match 
the traditional equivalent areas, which only include monopole (volume) and dipole (lift) effects. Analysis 
results using these traditional equivalent areas differ from the results obtained using off-body pressure 
distributions. An example is depicted in Figure 1, where the ground signatures calculated using the 
traditional equivalent area and the off-body pressure profile are superimposed and plotted. It is seen that 
for the same geometry, using the traditional Mach-cut equivalent area produces a ground signature that 
differs from that produced using off-body pressure distribution. This is especially significant considering 
that the accepted high-fidelity approach for sonic boom prediction on the ground involves propagation of 
the off-body pressure distribution. 

Even though traditional equivalent areas are attractive for their intuitive elegance and their reduced 
computational cost, given the above shortcoming in their analysis capability, they are not very useful in 
high fidelity shape optimization. To overcome this, a reversed equivalent area approach was developed 
and its design application demonstrated in earlier work. 1 ' The reversed equivalent area contains more 
information about the three dimensional flow for sonic boom analysis and produces essentially the same 
ground signature when a given geometry is analyzed by propagation using the off-body pressure distribu- 
tion. This is depicted in Figure 2, where the ground signatures propagated from the reversed equivalent 
area and the CFD off-body pressure distribution show excellent agreement. Using reversed equivalent 
area transforms equivalent area analysis to be on par with the off-body pressure distribution approach. 
Other researchers have used direct conversion of off-body pressure to equivalent area and used that in 
adjoint-based shape optimization. 18 From the perspective of getting a sonic boom ground signature, 
the directly converted equivalent area and the reversed equivalent area are both equivalent and superior 
to the traditional equivalent area. The key difference is that the reversed equivalent area offers better 
one-to-one correlation with the location and impact of the components of the aircraft concept under 
the assumption that 2-D propagation equations hold true closer to the aircraft. In this work, we limit 
ourselves to the reverse equivalent area, although directly converted equivalent areas can be used as well. 
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Figure 1. Shortcoming of Traditional Equivalent Areas 17 


This paper proposes using a discrete adjoint methodology to help the designer generate aircraft outer 
mold lines that attempt to achieve target reversed equivalent areas. Formal coupling of a CFD ad- 
joint methodology with a reverse boom propagation adjoint method, similar to that presented in earlier 
work, 16,9 allows efficient computation of the sensitivity of a reversed equivalent area based cost functional 
with respect to the aircraft shape design variables. The paper is organized as follows: Section II will 
provide a detailed mathematical derivation of the reversed equivalent area discrete adjoint formulation. 
Section III defines the problem setup including target equivalent area generation and surface parameter- 
ization. Section IV presents the optimization results obtained and refines the results with an alternative 
cost functional for boom mitigation. Section V will provide concluding remarks. 
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Figure 2. Advantage of Reversed Equivalent Areas 17 


II. Mathematics of the Reversed Boom Adjoint 

This section presents the mathematics behind the reversed boom adjoint methodology. The primal 
problem refers to the reversed augmented Burgers’ propagation, 17 fisted in Equation 1. 

dP _ dP ld 2 P C v gs 1 dG 1 d(poco) 

da dr T dr 2 "l + O^-^p 2G da 2p 0 co da 

An operator splitting scheme 19, 20 is used to solve a set of five equations under the assumption that 
if the time step is small, the error induced by splitting is small. As for the primal problem during 
propagation, the reversed propagation follows the same numerical steps. Equation 2 represents the effect 
of first relaxation and scaling due to ray tube area ( G ) spreading and stratification. The matrices included 
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in these equations are provided in the Appendix, and differ from the regular boom propagation problem 
mainly because of the change in propagation direction and the presence of regularization terms. The 
exact form and nature of these regularization terms is laid out in detail in earlier work. 1 ' Based on the 
discretization scheme used, the matrices are tridiagonal; hence the Thomas algorithm 21 is used to solve 
the system efficiently. Since there are two relaxation phenomena corresponding to oxygen and nitrogen, 
Equations 2 and 3 are each solved using their respective values for C v and 9 

A n q n = (2) 

A%r n = B?q n (3) 

For the absorption equation, a Crank-Nicholson scheme is used for advancing the pressure in time. 
Using this discretization scheme, the absorption phenomenon, when discretized, also transforms into a 
tridiagonal matrix problem as given in Equation 4, which is solved to obtain t n . 

A™t n = B™r n (4) 

The non-linear portion of Equation 1 is solved using an upwind flux-splitting method. This is written 
as shown in Equation 5, and the flux-splitting approximation is given in Equation 6, with Equations 7 
and 8 being the flux terms. 
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Expanding the terms results in the discretized equation for the non-linear part of the reversed Burgers’ 
equation primal problem as given in Equation 9. 
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The ray tube spreading and atmospheric stratification are simply scaling terms — these are included 
in the k factor in Equation 2. For the solution of the reversed augmented Burgers’ equation, Equations 2, 
3, 4 and 9 are solved repeatedly, in that order, for n = 1 . . . N time steps, and at each stage the pressure 
is updated, while also successively updating intermediate values r, q, and t. 


Adjoint 

The discrete adjoint equations for reverse propagation derived in this section are based on a similar 
implementation described in previous work. 9, 16 The adjoint equations are derived in two steps. In the 
first step, the sensitivity of the cost functional with respect to the off-body pressure is determined. In the 
second step, a formal coupling of reverse propagation and CFD allows computation of the sensitivity of 
the cost functional with respect to the aircraft shape parameters. In the following paragraphs, a detailed 
derivation of the first step is provided, followed by a brief description of the formal coupling process. 

A Lagrangian is first written to account for the reverse propagation process. If D is the vector of 
design variables (off-body dp/p in the first step) and I n is the cost functional, then the Lagrangian 
corresponding to this objective is given in Equation 10. Taking the derivative of the Lagrangian with 
respect to D results in Equation 11, where it is assumed that the cost functional does not depend explicitly 
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on the intermediate pressure vectors r, q , and t. Furthermore, the matrices themselves do not vary with 
the initial pressure profile. 
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Collecting the and terms from Equation 11 and equating them to zero results in 

four adjoint equations that are iteratively solved backwards in time. Collecting all the terms, and 
simplifying yields Equation 12. Similarly, collecting the and terms, we have Equations 

13, 14 and 15 respectively. The adjoint solution process involves solving Equations 12, 13, 14, and 15 
iteratively. Equation 12 is solved initially by assuming 7 o,iv+i = 0 since there are no N + 1 terms in our 
primal propagation problem. The intermediate adjoints are successively updated and solved. The primal 
problem is solved first, and relevant pressure vectors are stored for use in the adjoint process. 
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Previous studies have looked at specifying ground 16 as well as off-body t 8,22 targets to perform 
adjoint-based shape optimization. As mentioned in the introduction of this paper, equivalent areas offer 
a degree of intuition making them attractive to designers. The cost functional used in this formulation 
is the reversed equivalent area matching as given in Equation 16. The target reversed equivalent area 
is generated using spline and Bezier fits and is described in earlier work . 23 The derivative of the cost 
functional (Equation 17) can be used in Equation 12 to start the adjoint calculation process. However, 
the partial derivative of the reversed equivalent area with respect to the off-body pressure is needed. In 
order to get this term, conversion from pressure to E-function and from E-function to reversed equivalent 
area are considered. 
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According to Whitham , 24 the conversion from F-function to pressure is given by the simple scaling 
expression of Equation 18. The reversed equivalent area can be computed from the reverse propagated 
F-function values using Equation 19. When numerical integration is carried out, Equation 19 is recast as 
the summation equation given in Equation 20, with E& = OVfc < 0, Xk = 0 Vfc < 0. Using Equation 20, 
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the derivative term can be computed as given in Equation 21. Based on the definition of the reversed 
equivalent area in Equation 20, the derivative matrix a gp r is a square, lower-triangular matrix. Finally, 
Equation 18 is used to result in d gp r , which is then substituted in Equation 17. 
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Equation 9 is differentiated to obtain the partial derivative terms needed in the adjoint calculation. 
Taking the partial derivatives with respect to i”_ 1; f" and t r ‘ +1 yields Equations 22, 23 and 24 respectively. 
These are used to populate the Jacobian matrix in Equation 13. 
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After Equations 12-15 are solved iteratively, the gradient vector of the cost functional for the reverse 
propagation process is given by Equation 25. A formal process of coupling this with CFD is given in 
detail in earlier work. 9 Briefly stated, a set of three adjoint equations are solved for the interface between 
CFD and reverse propagation. The first equation uses the gradient vector computed in Equation 25 
to determine the boom interface [See Equation 22 in Reference 9 ] Lagrange multipliers. The other two 
adjoint equations solve for the Lagrange multipliers associated with CFD flow solution and mesh vectors 
respectively. Once solved, the sensitivity of the cost functional with respect to all the aircraft shape 
parameters is available for use by a gradient-based optimizer. 


— = Aar (25) 

The adjoint sensitivites obtained above are compared against those obtained using a complex step 
approach. The reverse propagation adjoint sensitivities as well as coupled-adjoint sensitivities match 
their complex step counterparts up to 13 digits of numerical precision. This verifies that the adjoint 
sensitivities are accurate in the numerical sense and can be used in the design process. 
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III. Problem Setup 


In this paper we optimize the baseline configuration shown in Figure 3. This configuration is the result 
of earlier optimization using mixed-fidelity, 25 non- adjoint-based reversed equivalent area 17 methods. The 
initial mesh for this concept was generated using VGRID, 26 SSGRID 27 and is shown in Figure 4. This grid 
generation approach is a heuristic technique for aligning the mesh topology a priori with the expected 
primary off-body shock structures. A more rigorous adjoint-based approach to mesh adaptation for such 
problems is described in literature. 11 The CFD grid utilizes a plane of symmetry along the centerline 
and contains 4 million nodes and 24 million tetrahedral elements. The surface mesh for the aircraft has 
been parameterized using a free-form shape deformation tool called BANDAIDS. 28 BANDAIDS provides 
a compact set of design variables for modifying a discrete surface mesh in the normal direction along with 
analytic sensitivities required by the discrete adjoint formulation of the near-field CFD problem. All the 
components of the baseline, except the flow through nacelle and pylon, are parameterized and allowed 
to vary in the shape optimization exercise. The intersections between aircraft components are held fixed 
for simplicity, although this is not a requirement of the formulation. A total of 138 design variables were 
used to parameterize the above-mentioned components of the aircraft concept, but only 81 of them are 
active during the optimization. 

L. 

Pod_X_S - 95.6 ft 

Pod_X_E~ 127.0 ft 


Fuselage_X_E - 117.0 ft 






Figure 3. Orthogonal projections of the baseline configuration 

The optimization problem is given in Equation 26. The reversed equivalent area corresponding to 
the baseline is obtained after computing the off-body flow solution. A combination of Bezier and spline 
curve fits are used to generate a smooth equivalent area profile that acts as an inverse- design target for 
the adjoint-based optimization. Because of the adverse sensitivity of deviations above the target and 
favorable sensitivity of deviations below the target, only perturbations above the target are penalized, 
while small perturbations below the target are accepted without counting towards the cost functional. 

1 M 

Minimize I — — £ [Aer.-A^r] (26) 

1=1 

IV. Optimization Results 

This section presents the shape optimization results of the concept in Figure 3 at a free-stream Mach 
numnber of 1.6 and angle of attack of 0.6 degrees. As viscous effects are likely to be small for these 
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Figure 4. CFD stretched grid 


configurations when the primary objective is sonic boom, an Euler solver is used for this study. The 
calculations were carried out on 22 nodes of the Altix ICE X cluster at the NASA Langley Research 
Center. Each node contains a dual socket, 8-core 2.6 Ghz Intel Sandy Bridge chipset amounting to 16 
cores per node and 32 GBs of memory (or 2 GBs per core). Different optimization packages such as 
NPSOL, 29 PORT, 30 and KSOPT 31 were tried during the course of this study. For the problem at hand, 
PORT was found to offer the best mix of performance while iterating towards the desired goal. Hence, 
PORT was used in the reversed equivalent area matching optimization presented in this study. 

Unlike off-body dp/p matching, which is localized in shape perturbations, shape changes for equiv- 
alent area matching tend to have a more global effect. As in off-body dp/p matching, equivalent area 
matching in the front portion is fairly straight-forward given that just the fuselage nose needs to be 
changed. However, when it comes to matching the aft portions of the equivalent area, all the components 
longitudinally ahead of the desired matching region, including the fuselage nose, will have an influence 
due to the hyperbolic nature of the flow equations in the supersonic regime. If all the shape parameters 
are allowed to vary at the same time, conflicting sensitivities of the cost functional with respect to certain 
shape parameters may cause the optimizer to make little progress or drive the optimizer away from the 
optimum in other regions. For example, if the optimizer is trying to match a target equivalent area 
that differs from the baseline distribution in both the front and aft sections, increasing the fuselage nose 
diameter at a particular section may get the distribution closer to the target in the front section and 
decreasing it may get it closer to the target in the aft sections. This conflict may cause the optimizer to 
stall before reaching an optimum. To overcome this, a multi-step optimization approach is used in this 
study to match the target equivalent area distribution. First, the nose is optimized to match the target 
distribution in the front portion. Then, the nose is frozen, and the mid-section is optimized. Finally 
the nose and mid-section are frozen and the aft is optimized. This forces the optimizer to reach the 
desired target reversed equivalent area distribution without destroying the matching obtained in previous 
optimizations in the march towards the aft regions. 

Figure 5 depicts the reversed equivalent area of the baseline and target overlaid with the reversed 
equivalent area after the multi-step adjoint-based shape optimization process. The adjoint-based shape 
optimization effectively shapes the baseline at appropriate regions to match the target area distribution. 
The final design equivalent area closely follows the target distribution till roughly 120 feet. Beyond 
this, the optimized concept equivavent area seems to oscillate slightly about the target equivalent area 
distribution. One or two more additional iterations that freeze the optimized components and march 
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towards shape optimization in the un-matched aft regions can perhaps improve the match. 



Figure 5. Comparison of Reversed Equivalent areas after Optimization 


Figure 6 shows the orthographic views of the baseline and equivalent-area matched designs. The 
parameterization chosen for this study allows smooth changes to the geometry with the maximum shape 
perturbation being 6 inches. Subtle changes to the nose, along with appreciable changes to the aft, espe- 
cially the pod (the non-lifting fuselage like component at the intersection of the vertical tail and horizontal 
tail), cause a change in the off-body pressure profile that generates a favorable reversed equivalent area 
distribution. 



Figure 6. Different views of the baseline and optimized concepts 


Figure 7 depicts the iteration history of the cost functional against the design cycle in each of the three 
optimizations carried out for equivalent area matching. For nose optimization, the optimizer ran 31 flow 
and 4 adjoint solutions. The optimization is terminated when the cost functional does not change in four 
consecutive design cycles. For the mid-section, the optimization gradually reduces the cost functional and 
terminates after 70 flow and 33 adjoint solutions. In the aft section, the cost functional drops significantly 
in the first few design cycles followed by slow and gradual reduction as the optimization progresses. For 
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this final phase, the optimization uses 33 flow and 23 adjoint solutions. For each of these cases, the 
flow solver was allowed to run for 800 iterations, during which the residual dropped by at least 8 orders 
of magnitude. The adjoint solver residual dropped typically by twelve orders of magnitude. On the 
hardware described above, each flow solve takes 4.5 minutes and each adjoint solution takes 3 minutes. 
In total, the optimization consumed roughly 13 hours of wall-clock time. 
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Figure 7. Iteration history of the optimization steps 

After shape optimization, a visual check of the perturbed geometries is performed to make sure that 
the optimizer does not generate any unreasonable cross-sections. Figures 8 shows the fuselage cross- 
sections, blown-up to show clarity, at the non-dimensional longitudinal stations indicated. The following 
observations can be made about the changes to the fuselage cross-sections: 

• The intersections are not allowed to vary. This is apparent from the cross-section at the longitudinal 
non-dimensional distance of ij = 0.8 

• The shape changes are smooth 

• Even though the baseline has circular cross-sections, the final design can have non-circular cross- 
sections, 


The aft fuselage is shrunk radially. 




Figure 8. Comparison of fuselage cross-sections 


Figures 9 and 10 show the comparison of the sections of the wing and horizontal tail respectively. 
For the wing, the thickness is fairly well maintained with very little changes to the top surface. Going 
from the in-board to out-board sections, there is a low-frequency undulation in the bottom surface that 
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disrupts a strong under-wing shock into multiple small shocks in the near-field. For the horizontal tail, 
the thickness is reduced for the out-board sections as the optimizer tries to match the target equivalent 
area distribution. The overall changes are smooth and free of artifacts. 



Figure 9. Comparison of wing cross-sections 




Figure 10. Comparison of horizontal tail cross-sections 


Changes to the pod geometry, depicted in Figure 11, produce significant changes to the equivalent 
area in the aft. The optimizer pinches the pod near the front while creating an aft bump. These changes, 
along with interaction with other components in the aft, take the design much closer to the desired target 
distribution. Figure 12 shows the changes to the vertical tail sections. The leading and trailing edge are 
not allowed to vary during optimization as seen from the sectional changes. The optimizer reduces the 
thickness of the vertical tail to compensate for the volume and lift changes created by changes to the 
pod, horizontal tail and the aft fuselage. 

Figure 13 shows the comparison of the ground signatures. The baseline has a perceived level of 80 
PLdB. The target equivalent area produces a smooth ground signature with a perceived level of 66.3 
PLdB. After the adjoint-based reversed equivalent area matching, the ground signature corresponding to 
the final design has a perceived level of 75.9 PLdB. Most significantly, the mid-shock is eliminated, while 
the front and aft portion shaping is improved. The small shocks seen in the ground signature are the 
result of equivalent area oscillations of the final design around the target equivalent area. 

Figure 14 shows the comparison of the baseline and final off-body pressure waveforms. The fuselage 
nose changes produce an oscillatory behavior that produces a smoother ground signature. The two-shocks 
generated by the wing of the baseline concept are broken into multiple smaller strength shocks in the 
final design. The aft shock system strengths and locations are subtly modified by the changes in the aft 
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Figure 11. Comparison of pod cross-sections 


components such that the ground signature is significantly better shaped in the aft. While none of the 
shocks in the aft are eliminated, their strengths are greatly reduced and their locations favorably altered. 



Figure 12. Comparison of4 vertical tail cross-sections 

The final design from the equivalent area matching was used as the starting point to minimize the 
ground A-weighted loudness cost functional. 9 Figure 15 shows that the front part of the ground signature 
is further smoothed by this process. This reduces the perceived level of loudness to 75.5 PLdB and the A- 
weighted loudness to 61.3 dBA. For the purpose of demonstrating the effectiveness of adjoint-based design 
optimization for matching reversed equivalent areas, additional iterations were not deemed necessary. It 
is the author’s belief that the approach presented in this paper perfectly complements other existing 
adjoint-based optimization formulations for sonic boom mitigation. Even though the demonstration 
case in this study was restricted to unconstrained under-track sonic boom mitigation of a flow-through 
concept, the approach can be a practical and useful tool for constrained under-track and off-track sonic 
boom optimization with powered engine boundary conditions. 

V. Conclusions 

The reversed equivalent area matching approach presented in this paper offers a third approach 
for carrying out adjoint-based shape optimization for sonic boom mitigation; the other two approaches 
use near-field and ground-based cost functionals respectively. Each of these approaches looks at the cost 
functional from a different perspective, and offers an efficient way to use adjoint-based shape optimization 
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Figure 14. Comparison of off-body pressure waveforms 
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as the concept proceeds from conceptual design to preliminary design. In the conceptual stages, the 
equivalent area matching presented in this paper may be used to not only deform the outer mold line, but 
also to locate components such as nacelles and control surfaces without adversely affecting the quality 
of computed gradients to reach a satisfactory level of matching. Equivalent area targets also provide 
guidance to the optimizer about the volume and lift changes needed to generate a low boom concept 
while allowing for constraints such as cockpit and cabin volume to be imposed. Finally, the ground 
loudness or near- field adjoint-based shape optimization may be used to further refine the concept. 



Figure 15. Comparison of ground signatures after loudness optimization 
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Appendix 

The tridiagonal matrices for the relaxation processes are: 
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B n ,B? 
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In the above matrices, Ki = — | C ^a" , «2 = = ~ AjC ’ = 2ii and a' = 1 — a. If using 

the Crank-Nicholson scheme, a = 0.5. For thermo- viscous absorption, the matrices are given below with 
A = — 2r ^") 2 > ^3 = ~ 2 aj 2 an d = 2^3- Ifi and f? 2 are the regularization parameters for the relaxation 
and absorption equations respectively. 


A% = 


( 1 0 ••• \ 
—A — t 3 (1 + 2A + £4) —A — t% • • • 


V 


Bgm 


( 1 0 

A + 1 3 (1 — 2A — £4) A + £3 


V 


0 1 

\ 


0 1 / 


16 of 16 


American Institute of Aeronautics and Astronautics 


